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Abstract. A number of astrophysical scenarios possess and preserve an overall cylindrical 
symmetry also when undergoing a catastrophic and nonlinear evolution. Exploiting such 
a symmetry, these processes can be studied through numerical-relativity simulations at 
smaller computational costs and at considerably larger spatial resolutions. We here present 
a new flux-conservative formulation of the relativistic hydrodynamics equations in cylindrical 
coordinates. By rearranging those terms in the equations which are the sources of the largest 
numerical errors, the new formulation yields a global truncation en'or which is one or more 
orders of magnitude smaller than those of alternative and commonly used formulations. We 
illustrate this through a series of numerical tests involving the evolution of oscillating spherical 
and rotating stars, as well as shock- tube tests. 
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1. Introduction 

Numerical simulations assuming and enforcing axisymmetry are particularly useful to study 
at higher resolution and smaller computational costs those astrophysical scenarios whose 
evolution is expected to possess and preserve such a symmetry. On the other side, the 
numerical solution of systems of equations expressed in coordinates adapted to the symmetry 
has often posed serious difficulties, because of the coordinate singularity present on the 
symmetry axis. The "cartoon" method, proposed by Alcubierre et al. [1], allows to exploit 
the advantages of reduced computational resource requirements, while adopting Cartesian 
coordinates, which are non-singular 

The "cartoon" method proves particularly useful in the numerical evolution of smooth 
functions, like the metric quantities of the Einstein equations. However, because of the 
interpolations necessary to impose the axisymmetric conditions on a Cartesian grid, the 
"cartoon " approach is not considered to be accurate enough to describe the shocks which 
generically develop when matter is present. As a consequence, general-relativistic codes 
employing the "cartoon" method have adopted cylindrical coordinates for the evolution of 
the matter (and magnetic field) variables [2, 3, 4, 5, 6]. All the cited works adopt the 
same formulation for the hydrodynamical equations in cylindrical coordinates. In the present 
article, we propose a slightly different formulation, which has proven to reduce the numerical 
errors, especially in the vicinity of the symmetry axis. 

More specifically, we have written the Whisky2D code, which solves the general-relativistic 
hydrodynamics equations in a flux-conservative form and in cylindrical coordinates. This of 
course brings in l/r singular terms, which must be dealt with appropriately. In the above- 
referenced works, the flux operator is expanded and the l/r terms, not containing derivatives, 
are moved to the right hand side of the equation (the source term), so that the left hand side 
assumes a form identical to the one of the three-dimensional (3D) Cartesian formulation. 
We call this the standard formulation. An other possibility is not to split the flux operator 
and to redefine the conserved variables, via a multiplication by r. We call this the new 
formulation. The new equations are solved with the same methods as in the Cartesian case. 
From a mathematical point of view, one would not expect differences between the two ways 
of writing the differential operator, but, of course, a difference is present at the numerical 
level. Our tests show that the new formulation yields results with a global truncation error 
which is one or more orders of magnitude smaller than those of alternative and commonly 
used formulations. 

Here we perform a series of tests to ascertain the convergence behaviour of the two 
formulations. We then show that the new formulation produces results which are generally 
more accurate, with a truncation error which can be several orders of magnitude smaller 
The paper in organized as follows. In Section 2 we remind the essentials of the "cartoon" 
approach for the evolution of the geometrical variables, while in Section 3 we review the 
flux-conservative formulation of relativistic hydrodynamics. We write down the relativistic 
flux-conservative hydrodynamics equations for axisymmetric formulations and we illustrate 
the two possible ways to write the singular term. In Section 4, we present several tests 
that compare the two formulations. We begin with the conservation of rest mass and 
angular momentum in the Cowling approximation and in full-spacetime evolution. Then the 
eigenfrequencies of uniformly rotating neutron-star models are compared with the results of 
a perturbative code. The last test examines the differences between the two formulations with 
respect to an analytic solution of an extreme shock case, which mimics the reflection of a cold 
and very fast gas at the symmetry axis. 

We have used a spacelike signature (— , +, +, +), with Greek indices running from to 3, 
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Latin indices from 1 to 3 and the standard convention for the summation over repeated indices. 
Unless expUcitly stated, all the quantities are expressed in the system of dimensionless units 
in which c = G = Mq = 1 . 



2. Evolution of the Einstein equations 

The logical and algorithmic structures of the Whisky2D code presented here follow closely 
the ones of the CCATIE [7] and Whisky codes [8], which solve the same set of equations in 
3D and using Cartesian coordinates. In what follows we provide only a brief overview of the 
set of equations for the evolution of the fields (within the "cartoon" prescription [1]) and for 
the evolution of the fluid variables, referring the interested reader to refs. [7, 9, 10] for a more 
detailed discussion. As for the other codes mentioned above, also Whisky2D is based on the 
Cactus Computational Toolkit [11]. 

More specifically, we evolve a conformal-traceless "3 + 1" formulation of the Einstein 
equations [12, 13, 14], in which the spacetime is decomposed into 3D spacelike slices, 
described by a metric jij, its embedding in the full spacetime, specified by the extrinsic 
curvature Kij, and the gauge functions a (lapse) and (shift), which specify a coordinate 
frame (see [15] for a general description of the 3 + 1 split). The particular system which we 
evolve transforms the standard ADM variables as follows. The 3-metric -fij is conformally 
transformed via 

1 



12 



In det 7ij , 



-4* 



7y, 



(1) 



and the conformal factor $ evolved as an independent variable, whereas 7^ is subject to 
the constraint dct7y = 1. The extrinsic curvature is subjected to the same conformal 
transformation and its trace tr Kij is evolved as an independent variable. That is, in place 
of Kij we evolve: 



K = tr K, 



= .9'^ Ay, 
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with tr A. 



0. Finally, new evolution variables 
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are introduced, defined in terms of the Christoffel symbols of the conformal 3-metric. 

The Einstein equations specify a well known set of evolution equations for the listed variables 

and are given by 
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where i?y is the three-dimensional Ricci tensor, Di is the covariant derivative associated with 
the three metric 7^ , "TF" indicates the trace-free part of tensor objects and Padm' ^^'^ ^ij 
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are the matter source terms defined as 

Padm = ria^npT"''^, 

S, = --f^c.n|3T"'^ , (9) 

where Ua = (—a, 0, 0, 0) and T"'^ is the stress-energy tensor for a perfect fluid (see 
Section 3). 

Four elliptic constraint equations, which are usually referred to as Hamiltonian and 
momentum constraints, 

n = +K^~ K,jK'' - 16^p^„^, = , (10) 
= Dj{K'^ - f^K) - 8ttS' = , (11) 

should be satisfied within each spacelike slice. Here R'-^^ = Rij'j^^ is the Ricci scalar on a 
3D timeslice. Additional constraints are given by 

det7,, = l, trl,, -0, r = f^f)^ ^ (12) 

with the last two equations of (12) being enforced algebraically. The remaining constraint 
in (12) and the constraints H and A^* are not actively enforced and can be used as monitors 
of the accuracy of our numerical solution. 

We specify the gauges in terms of the standard ADM lapse function, a, and shift vector, 
/3" [16]. We evolve the lapse according to the "1 + log" slicing condition [17]: 

dta - P'd.a = -2a{K - Kq), (13) 

where Kq is the initial value of the trace of the extrinsic curvature and equals zero for 
the maximally sliced initial data we consider here. The code uses a hyperbohc F-driver 
condition [18] 

dtf3'-f3^dj(3' =laB\ (14) 

dtB'- (3^ djB' ^ dt r - 13^ djr~r]B\ (15) 

where ij is a parameter which acts as a damping coefficient (see discussion in ref. [19]). 
Two routes are possible when solving numerically the Einstein equations in axisymmetric 
spacetimes. One route consists in using coordinates that exploit the symmetry and enforce its 
preservation already at a mathematical level, such as cylindrical coordinates. This advantage 
is counterbalanced by the fact that such coordinates are usually singular somewhere (e.g., on 
the axis for cylindrical coordinates) and that regularization conditions are therefore necessary 
(see [20, 21] and references therein for a recent discussion). 

The second route consists, instead, in using Cartesian coordinates and in exploiting the 
fact that these coincide with the cylindrical ones in one plane, namely the {x, z) plane 
(for concreteness we will assume hereafter that the Cartesian and the cylindrical z-axes 
coincide). The chief advantages of this approach, which is usually referred to as the "cartoon " 
method [1], are the absence of the need of regularization conditions and the easiness of 
implementation, through a simple dimensional reduction from fully 3D codes in Cartesian 
coordinates. However, these advantages are counterbalanced by at least two disadvantages. 
The first one is that the method still essentially requires the use of a 3D domain covered with 
Cartesian coordinates, although one of the three dimensions, namely the y-direction, has a 
very small extent. The second one is that, in order to compute the spatial derivatives in the 
2/-direction appearing in the Einstein equations, a number of high-order interpolations onto 
the X-axis are necessary (see discussion below) and these can amount to a significant portion 
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of the time spent for each evolution to the new timelevel. In practice, the spatial derivatives 
in the y-direction are computed exploiting the fact that all quantities are constant on cylinders 
and thus the value of a variable at a generic position (x, y, z) off the (x, z) plane can be 
computed from the corresponding value ^'(i, 0, z) on the (a;, z) plane, where 

x = + (^g) 

Clearly, since the solution of the evolution equations is computed only on the {x,z) plane, 
interpolations (with truncation errors smaller than that of the finite-difference operators) are 
needed at all the positions (i, y = 0, z). 

Overall the "cartoon" method represents the choice for many codes and it has been 
implemented with success in many applications, e.g., [1, 2, 3, 4, 22, 23, 5, 6] to cite a few. 



3. Evolution of the relativistic hydrodynamics equations 

An important feature of multidimensional non-vacuum numerical-relativity codes that solve 
the coupled Einstein-hydrodynamics equations in Cartesian coordinates is the adoption of a 
conservative formulation of the hydrodynamics equations [24, 25]. In such a formulation, 
the set of conservation equations for the stress-energy tensor T^'' and for the matter current 
density J'', that is 

V,,J^ = 0, (17) 

V^T^'' = 0, (18) 

is written in a hyperbolic, first-order "flux-conservative" form of the type [26] 

9tq + ftf«(q)=s(q), (19) 

where f (q) and s(q) are the flux vectors and source terms, respectively [27]. Note that the 
right-hand side (the source terms) depends only on the metric, on its first derivatives and on 
the stress-energy tensor Furthermore, while the system (19) is not strictly hyperbolic, strong 
hyperbolicity is recovered in a flat spacetime, where s(q) = 0. 

As shown by [25], in order to write the system (17)-(18) in the form of system (19), the 
primitive hydrodynamical variables {i.e. the rest-mass density p, the pressure p measured in 
the rest-frame of the fluid, the fluid 3-velocity measured by a local zero-angular momentum 
observer, the specific internal energy e and the Lorentz factor W) are mapped to the so called 
conserved variables q = {D, S"*, r) via the relations 

D=^pW, 

S' = y/jphW'^v' , (20) 
T =^{phW^~p)-D, 

where h = l + e + p/pis the specific enthalpy and W = (1 ~ 7y 

The advantage of a flux-conservative formulation is that it allows to use high-resolution shock- 
capturing (HRSC) schemes, which are based on Riemann solvers and which are essential for 
a correct representation of shocks. This is particularly important in astrophysical simulations, 
where large shocks are expected. In this approach, all variables q are represented on the 
numerical grid by cell-integral averages. The function is then reconstructed within each 
cell, usually through piecewise polynomials, in a way that preserves the conservation of the 
variables q. This gives two values at each cell boundary, which are then used as initial data for 
the (approximate) Riemann problem, whose solution gives the flux through the cell boundary. 
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As in the Whisky code, the evolution equations are here integrated in time using the 
method of Une [28], which reduces the partial differential equations (19) to a set of ordinary 
differential equations that can be evolved using standard numerical methods, such as Runge- 
Kutta or the iterative Cranck-Nicolson schemes [29, 30]. Furthermore, the Whisky2D 
code implements several reconstruction methods, such as Total- Variation-Diminishing (TVD) 
methods, Essentially-Non-Oscillatory (ENO) methods [31] and the Piecewise-Parabolic- 
Method (PPM) [32]. Also, a variety of approximate Riemann solvers can be used, starting 
from the Harten-Lax-van Leer-Einfeldt (HLLE) solver [33], over to the Roe solver [34] and 
the Marquina flux formula [35] (see [8, 9] for a more detailed discussion). 
The ability of properly evolving large gradients moving at relativistic speeds represents one 
of the main motivations that make this formulation the choice for all of the present 3D 
numerical-relativity codes solving the relativistic hydrodynamics equations on Eulerian grids 
(see refs. [36, 37, 38] for some of the most recent examples and ref. [39] for an alternative 
Lagrangian method). However, when considered within the axisymmetric approach used 
here, the use of a flux-conservative formulation in Cartesian coordinates may suffer from a 
potentially very serious disadvantage. In fact, the interpolations required by the "cartoon " 
method may be highly inaccurate when discontinuities in the fluid variables appear. To 
confront this problem, Shibata [2] has made the useful suggestion of writing the relativistic 
hydrodynamics equations in cylindrical coordinates, while keeping the solution of the Einstein 
equations in Cartesian coordinates. This approach has the obvious advantage that it does not 
require interpolation and that it exploits, at the mathematical level, the symmetries of the 
system, thus guaranteeing a better conservation of mass and angular momentum. However, 
the use of cylindrical coordinates for the evolution of the fluid variables also comes with an 
undesirable property; the coordinates are degenerate at the symmetry axis and the equations 
are no longer free of singularities. As we will comment in the following Section, this 
drawback can be compensated through a suitable formulation of the equations and a proper 
setup of the numerical grid. 



3.1. A new formulation of the hydrodynamics equations 

As mentioned in the previous Section, following ref. [3], we write the relativistic 
hydrodynamics equations (17)-(18) in a first-order form in space and time using cylindrical 
coordinates (r, 0, z). However, as an important difference from the approach suggested in 
ref. [3], we do not introduce source terms that contain coordinate singularities. Rather, we 
re-define the conserved quantities in such a way to remove the singular terms, which are the 
largest source of truncation error, also when evaluated far from the axis. 
We illustrate our approach by using as a representative example the continuity equation. This 
is the simplest of the five hydrodynamical equations but already contains all the basic elements 
necessary to illustrate the new formulation. We start by using the definitions for the conserved 
variables (20) to write eq. (17) generically as 

dtiV^lpW) + y^pW {av' - Z?')] = , (21) 
which in cylindrical coordinates takes the form 

dA^p"^) + dr \^pW {av- - 13-)] + [^pW {av' - /3^)] = , (22) 

where ^7 is the determinant of the 3-metric in cylindrical coordinates and where we have 
enforced the condition of axisymmetry d,p = 0. Because any 0-constant plane in cylindrical 
coordinates can be mapped into the {x, z) plane in Cartesian coordinates, we consider 
equation (22) as expressed in Cartesian coordinates and restricted to the y = plane, i.e., 

dt{xD) + [xD (aw^ - /?^)] + [xD {av' - f3')] = , (23) 
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where we have exploited the fact that for any vector of components A* on this plane 
= j4^, A'f' = Ay and 7 ~ x^^, with 7 being the determinant of the 3-metric in 
Cartesian coordinates. Equation (23) represents the prototype of the formulation proposed 
here, which we will refer to hereafter as the "new" formulation to contrast it with the 
formulation adopted so far, e.g., in ref. [3], for the solution of the relativistic hydrodynamics 
equations in axisymmetry and in Cartesian coordinates. The only, but important, difference 
with respect to the "standard" formulation is that in the latter the derivative in the x-direction 
is written out explicitly and becomes part of the source term s(q), i.e., 

D (aw^ - /?^) 



dt{D) + dx [D {av- - r )] + [D {av' - p')] 



(24) 



Even though the right-hand-side of eq. (24) is never evaluated at x = (because no grid 
points are located at a; = 0), both the numerator and the denominator of the right-hand-side of 
eq. (24) are very small for x ~ 0, so that small round-off errors in the evaluation of the right- 
hand-side can increase the overall truncation error. Stated differently, the right-hand-side of 
eq. (24) becomes stiff for x ~ and this opens the door to the problems encountered in the 
numerical solution of hyperbolic equations with stiff source terms [40]. 
What was done for the continuity equation (23) can be extended to the other hydrodynamics 
equations which, for the conservation of momentum in the x- and z-directions, take the form 



1 



dt {xSa) + dx [x {Sa [av- - ^ ) + a^p6\)] 



dAx{SA{av' -n + aV^pSl)] 



adAa ) + T'^'P'dAla + r°9^/3^ + -T'^'dAli. 



(25) 



with A = x,z. Similarly, the evolution of the conserved angular momentum S^ 
expressed as 

1 



ax^ 



dt (x^Sy) + dx [x^Sy {av- ~ r )] + d, [x^Sy {uv' - P')] \ = 



while the equation of the energy conservation is given by 
1 



ax^/j 



dt {xt) + dx [x (r {av- - ^ ) + pv-)] + d, [x (r {av' pv')] 



cS*,, is 



(26) 



(27) 



The changes made to the formulation are rather simple but, as we will show in Section 4, 
these can produce significant improvements on the overall accuracy of the simulations with a 
truncation error at least one order of magnitude smaller for all of the tests considered. Because 
of its simplicity, the changes in the new formulation of the equations can be implemented 
straightforwardly in codes written using the standard formulation. 

Finally, we note that both eq. (24) and eq. (23) are written in a flux-conservative form in the 
sense that the source term does not contain first-order spatial derivatives of the conserved 
variables. More precisely, eq. (23) is written in a flux-conservative form, while eq. (24) is 
written in a "flux-balanced" form, as it is typical for flux-conservative equations written in 
curvilinear coordinates [26] . The same is true also for eqs. (25)-(27) and for the corresponding 
equations presented in ref. [3], which are incorrectly classified as non flux-conservative. 
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Table 1. Equilibrium properties of the initial stellar models. The different columns refer 
respectively to: the ratio of the polar to equatorial coordinate radii Vp/re, the central rest-mass 
density pc, the gravitational mass _A/, the rest mass A/q, the circumferential equatorial radius 
Re, the angular velocity Q, the maximum angular velocity for a star of the same rest mass 
the ratio J/_A/^ where J is the angular momentum, the ratio of rotational kinetic energy 
to gravitational binding energy T/|iy|. All models have been computed with a polytropic 
EOS with K = 100 and T = 2. 
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Mo 


Re 
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T/\W\ 
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(Me) 
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A 


1.00 


1.28 


1.400 


1.506 


9.586 


0.000 


3.987 


0.000 


0.000 


B 


0.67 


1.28 


1.651 


1.786 


12.042 


0.253 


3.108 


0.594 


0.081 



3.2. Equation of state 

In whatever coordinate system they are written, the system of hydrodynamics equations can 
be closed only after specifying an additional equation, the equation of state (EOS), which 
relates the pressure to the rest-mass density and to the energy density. The code has been 
written to use any EOS, but all the tests so far have been performed using either an (isentropic) 
polytropic EOS 

P = Kp^, (28) 
e = P+ ^ , (29) 
or an "ideal-fluid" EOS 

p=iT-l)pe. (30) 

Here, e is the energy density in the rest frame of the fluid, K the polytropic constant (not 
to be confused with the trace of the extrinsic curvature defined earlier) and F the adiabatic 
exponent. In the case of the polytropic EOS (28), F = 1 + 1/N, where is the polytropic 
index and the evolution equation for r does not need to be solved. In the case of the ideal- 
fluid EOS (30), on the other hand, non-isentropic changes can take place in the fluid and the 
evolution equation for r needs to be solved. Note that the polytropic EOS (28) is isentropic 
and thus does not allow for the formation of physical shocks, in which entropy (and internal 
energy) can be increased locally (shock heating). 



4. Numerical tests 

In order to test the stability properties of the new formulation and compare its accuracy with 
that of the formulation first presented in [3] and then used, among others, in [4, 6, 41, 42], we 
have implemented both of them in Whisky2D. After this paper was published we were made 
aware via a private communication [43] that the numerical approach followed in [23] is in 
practice very similar to our new formulation, although the version in cylindrical coordinates 
of the equations discussed in [23] is not in a flux-conservative form; information about such 
numerical implementation was not given in [23] and therefore not available to us at the time 
this work was written. 

The initial data, in particular, has been produced as solution of the Einstein equations for 
axisymmetric and stationary stellar configurations [44], using the EOS (28) with F = 2 and 
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polytropic constant K — 100, in order to produce stellar models that are, at least qualitatively, 
representative of what is expected from observations of neutron stars. Our attention has been 
restricted to two illustrative models representing a nonrotating star and a rapidly rotating star 
having equatorial and polar (coordinate) radii in a ratio rp/re = 0.67. The relevant properties 
of these stellar models are reported in Table 1 . 

All the numerical results presented hereafter have been obtained with the following fiducial 
numerical set-up: the reconstruction of the values at the boundaries of the computational cells 
is made using the PPM method [32], while the HLLE algorithm is used as an approximate 
Riemann solver [33]. The lapse function is evolved with the "1 + log" slicing condition 
given by eq. (13), while the shift is evolved using a version of the hyperbolic F-driver 
condition (14) in which the advection terms for the variables and are set to zero. 

The time evolution is made with a method-of-line approach [28] and a third-order Runge- 
Kutta integration scheme (our CFL factor is usually chosen between 0.3 and 0.5). A third- 
order Lagrangian interpolation is adopted to implement the "cartoon " method. For the matter 
variables we use "Dirichlefboundary conditions (i.e., the solution at the outer boundary is 
always kept to be the initial one), while for the field variables we adopt outgoing Sommerfeld 
boundary conditions. 

We typically present results at four different resolutions: h — 0.4A/, h/2, h/i, 3h/16 and 
h/8, which correspond to about 25, 50, 100, 133 and 200 points across the stellar radius, 
respectively. The computational domain extends to 20 M both in the x and z directions, and 
a reflection symmetry is applied across the equatorial {i.e., z = 0) plane. Finally, we remark 
that in contrast with the interesting analysis of [45], we could not find signs of numerical 
instabilities when using the above numerical prescriptions for either of the two formulations 
considered. 

4.1. Oscillating Neutron stars: fixed spacetime 

The first set of tests we discuss has been carried out by simulating relativistic polytropic stars 
in equilibrium and in a fixed spacetime {i.e. in the Cowling approximation). In this case the 
Einstein equations are not evolved and the truncation error is in general smaller because it is 
produced uniquely from the evolution of the hydrodynamics equations. 
Although the stars are in equilibrium, oscillations are triggered by the first-order truncation 
error at the center and the surface of the star (our hydrodynamical evolution schemes are 
only first order at local extrema). Both the amplitude of the oscillations and the rate of the 
secular change in their amplitude converge to zero at nearly second order with increasing 
grid resolution [46, 47]. The genuine dynamics produced by the truncation error can then be 
studied either when the spacetime is held fixed {i.e., in the Cowling approximation) or when 
the spacetime is evolved through the solution of the Einstein equations. This is shown in 
Fig. 1, which reports the evolution of the central rest-mass density for rapidly rotating stars 
(model B in Table 1) evolved within the Cowling approximation. The left panel refers to the 
standard formulation, while the right one to the new formulation. Note that in both cases the 
amplitude of the oscillations decreases at roughly second order with increasing resolution, 
while keeping the same phase. This is a clear signature that the oscillations corresponds to 
proper eigenmodes of the simulated star. However, the difference of the secular evolution 
between the standard formulation and the new one is rather remarkable. The latter, in fact, is 
much more accurate and the well-known secular increase in the central density is essentially 
absent in the new formulation. 

Quantities that are particularly useful to assess the accuracy of the two formulations are the 
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Figure 1. Evolution of the central rest-mass density for rapidly rotating stars (model B in 
Table 1) evolved within the Cowling approximation. The left panel refers to the use of the 
standard formulation, while the right one to the new formulation. Note the different scales 
in the two panels and note that in both cases the amphtude of the oscillations decreases with 
increasing resolution, while keeping the same phase. 



rest mass and the angular momentum which we compute as [48] 



J. = 



2"" / y^pWx dx dz , 



(31) 
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xdxdz , (32) 



where K is the coordinate volume occupied by the star and V is coordinate volume of the 
computational domain. 

Figure 2 shows the dependence on the inverse of the resolution of the error in the conservation 
of the rest mass for a nonrotating model as computed in the Cowling approximation (left 
panel) or in a fully dynamical simulation (right panel). Since the evolution of the rest mass 
shows, in addition to a secular evolution, small oscillations (i.e., of ~ 3 x 10~^ for the 
highest resolution and of 3 x 10^® for the lowest resolution) the calculation of the rest 
mass at a given time can be somewhat ambiguous. To tackle this problem and to avoid 
the measurement to be spoiled by the oscillations, we perform a linear fit of the evolution 
of A/o, normalized to the initial value Mo{t = 0), between the initial value and a time 
t ~ 25 ms (corresponding to about 30 oscillations) and we take as the time derivative of 
the mass the coefficient of the linear fit: d{Mo / Mo{t = 0))/dt. Fig. 2, in particular, reports 
in a logarithmic scale d(A/o/.A/o(^ = 0))/d^ a function of the inverse of the resolution 
h. Indicated with squares are the numerical values obtained with the standard formulation of 
the hydrodynamics equations, while triangles are used for the new one. Also indicated with a 
long-short-dashed line is the slope for a second-order convergence rate. 
Note that although we use a third-order method for the reconstruction (namely, PPM), we do 
not expect third-order convergence. This is also due to the fact that the reconstruction schemes 
are only first-order accurate at local extrema {i.e. at the centre and at the surface of the star), 
thus increasing the overall truncation error. Similar estimates were obtained also using the 
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Figure 2. Time derivative of the average of tlie rest mass Mo , normalized to tlie initial value 
Mo{t = 0), for evolutions in a fixed spacetime (Cowling approximation). The average 
d(Mo/A/o(t = 0))/dt is computed between the initial value and a time t = 25ms, 
coiTesponding to about 30 oscillations. The left panel refers to a nonrotating star (model A in 
Table 1), while the right panel to a rapidly rotating star (model B in Table 1). Indicated with 
squares are the numerical values obtained with the standard formulation of the hydrodynamics 
equations, while triangles are used for the new one. Also indicated with a dot-dashed line is 
the slope for a second-order convergence rate. 
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Figure 3. Time derivative of average of the angular momentum normalized to the initial value 
d( J/ J{t = 0))/dt (cf.. Fig. 2) for a rapidly rotating star (model B of Table 1). Indicated with 
squares are the numerical values obtained with the standard formulation of the hydrodynamics 
equations, while triangles are used for the new one; a dot-dashed line is the slope for a second- 
order convergence rate. 



Whisky code in 3D Cartesian coordinates [8, 9]. 

Clearly both the new and the standard methods provide a convergence rate which is close to 
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two. However, and this is the most important result of this work, the new method yields 
a truncation error which is several orders of magnitude smaller than the old one. More 
specifically, in the case of the rest mass, the conservation is more accurate of about four 
orders of magnitude. We believe that this is essentially due to the rewriting of the source 
terms in the flux-conservative formulation which in the new formulation does not have any 
coordinate-singular term (i.e. oc l/x). 

Note also that, because the new formulation is intrinsically more accurate, it also suffers more 
easily from the contamination of errors which are not directly related to the finite-difference 
operators. [The one made in the calculation of the integral (32) is a relevant example but it 
is not the only one]. This may be the reason why, in general, at lower resolutions the new 
formulation has convergence rate which is not exactly two and appears over-convergent (see 
right panel of Fig. 2). However, as the resolution is increased and the finite-difference errors 
become the dominant ones, a clearer trend in the convergence rate is recovered. 
Another way of measuring the accuracy of the two formulations is via the comparison of the 
evolution of the angular momentum. While this quantity is conserved to machine precision 
in the case of a nonrotating star, this does not happen for rotating stars and the error can be 
of a few percent in the case of very low resolution and of a very rapidly rotating star This is 
shown in Fig. 3 for the stellar model B of Table 1 and it reports in a logarithmic scale the time 
derivative of the average of the angular momentum J normalized to the initial value J{t = 0). 
In analogy with Fig. 2, in order to remove the small-scale oscillations we first perform a linear 
fit of the evolution of J between the initial value and a time t ^ 25 ms and take the coefficient 
of the fit as the time derivative of the angular momentum: d( J/ J{t = 0)) /dt. Indicated with 
squares are the numerical values obtained with the standard formulation of the hydrodynamics 
equations, while triangles are used for the new one; a dot-dashed line shows the slope for a 
second-order convergence rate. 

It is simple to recognize from Fig. 3 that also for the angular momentum conservation the new 
formulation yields a truncation error which is two or more orders of magnitude smaller, with 
a clear second-order convergence being recovered at sufficiently high resolution. 




Figure 4. The same as in Fig. 1 but for a full-spacetime evolution. Tlie left panel refers to 
the standard formulation, while the right one to the new formulation. Note the different scale 
between the two panels. 
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4.2. Oscillating Neutron stars: dynamical spacetime 

Also the second set of tests we discuss is based on the evolution of relativistic polytropic 
stars in equilibrium, but now the evolution is performed in a dynamical spacetime, thus with 
the coupling of Einstein and hydrodynamics equations. The truncation error in this case is 
given by the truncation error coming from the solution of both the field equations and the 
hydrodynamics equations. The results of our calculations are summarized in Figs. 4-6, which 
represent the equivalents of Figs. 1-3 for full-spacetime evolutions. Because the results are 
self-explanatory and qualitatively similar to the ones discussed for the evolutions with fixed 
spacetimes, we will comment on them only briefly. 




2 4 6 8 10 20 2 4 6 8 10 20 

1/h (M^i) 1/h (Mgi) 



Figure 5. The same as in Fig. 2, but for full-spacetime evolutions. The left panel refers to a 
noni'otating star (model A in Table 1 ), while the right panel to a rapidly rotating stai' (model B 
in Table 1). 

In particular. Figs. 5-6 highlight that while the overall truncation error in dynamical 
spacetimes is essentially unchanged for the standard formulation, it has increased in the case 
of the new formulation. This is particularly evident at very low resolutions, where the new 
formulation seems to be hyper-convergent. However, despite a truncation error which is larger 
than the one for fixed spacetimes, the figures also indicate that the new formulation does 
represent a considerable improvement over the standard one and that its truncation error is 
at least two orders of magnitude smaller. Most importantly, the conservation properties of 
the numerical scheme have greatly improved and the secular increase in the rest mass, is 
also considerably suppressed. This is clearly shown in Fig. 4, where the secular increase 
is suppressed almost quadratically with resolution. More precisely, for both approaches the 
growth rate of the central rest-mass density for the coarse resolution is ~ 12 times larger than 
the corresponding one for the high resolution. However, at the highest resolution, the growth 
rate for the standard formulation is ^ 10 times larger than the one of the new formulation. 

4.3. Calculation of the eigenfrequencies 

As mentioned in the previous Section, although in equilibrium, the simulated stars undergo 
oscillations which are triggered by the nonzero truncation error It is possible to consider 
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Figure 6. The same as Fig. 3 but for a rapidly rotating star evolved in a dynamical spacetime. 



these oscillations not as a numerical nuisance, on the contrary it is possible to exploit 
them to perform a check on the consistency of a full nonlinear evolution with a small 
perturbation (the truncation error) with the predictions of perturbation theory [46, 47]. 
Furthermore, when used in conjunction with highly accurate codes, these oscillations can 
provide important information on the stellar oscillations within regimes, such as those of very 
rapid or differential rotation, which are not yet accessible via perturbative calculations [49]. 
In this Section we use such oscillations, and in particular the fundamental £ = quasi-radial 
f -mode, to compare the accuracy of the two formulations against the perturbative predictions. 
This is summarized in Fig. 7 which reports the power spectral density (in arbitrary units) of the 
maximum rest-mass density evolution (c/. Figs. 1 and 3) in the new and standard formulation 
(solid and dashed lines, respectively). The simulations are relative to a nonrotating star (model 
A in Table 1) with the left panel referring to an evolution with a fixed spacetime, while 
the right one to an evolution with a dynamical spacetime. The specific spectra shown are 
calculated from the simulations at the highest resolution and cover an interval of 25 ms. It 
is quite apparent that the two formulations yield spectra which are extremely similar, with 
a prominent F-mode at about 2.7 kHz and 1.4 kHz for the fixed and dynamical spacetime 
evolutions, respectively. The spectra also show the expected quasi-radial overtones at roughly 
multiple integers of the F-mode, the first of which has a comparable power in the case of 
Cowling evolution, while it is reduced of about 50% in the full spacetime evolution. Indeed, 
the spectra in the two formulations are so similar that it is necessary to concentrate on the 
features of the F-mode to appreciate the small differences. These are shown in the insets of 
the two panels which report, besides a magnification of the spectra near the F-mode, also the 
perturbative estimate Fport, as calculated with the perturbative code described in ref. [50]. 
To provide a more quantitative assessment of the accuracy with which the two formulations 
reproduce the perturbative result we have computed the eigenfrequency of the F-mode, which 
we indicate as Fnum, by performing a Lorentzian fit to the power spectrum with a window of 
0.2 kHz. We remark that it is essential to make use of a Lorentzian function for the fit as this 
reflects the expected functional behaviour and increases the accuracy of the fit significantly. 
Shown in Fig. 8 is the absolute value of the relative difference between the numerical and 
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Figure 7. Power spectral density (in arbitrary units) of tlie maximum rest-mass density 
evolution in tlie new and standard formulation (solid and dashed lines, respectively). The 
simulations are relative to a nonrotating stai' (model A in Table 1 ) with the left panel referring 
to an evolution with a fixed spacetime and the right one to an evolution with a dynamical 
spacetime. The spectra are calculated from the simulations at the highest resolution and cover 
25 ms of evolution. For both panels the insets show a magnification of the spectra near the 
F-mode and the comparison with the perturbative estimate as calculated with the numerical 
code described in ref. [50]. 




Figure 8. Relative difference between the numerical and perturbative eigenfrequencies of the 
F-mode for the two formulations (solid lines for the new one and dashed lines for the standard 
one). The differences are computed for different resolutions and refer to the nonrotating model 
A of Table 1 when evolved in a fixed spacetime (left panel) and in a dynamical one (right 
panel). Indicated with a dot-dashed line is the slope for a second-order convergence rate. 



perturbative eigenfrequencies of the F-mode, 1 1 — Fnum/ ^pert | for the two formulations (solid 
Unes for the new one and dashed lines for the standard one). The differences are computed 
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for different resolutions with h = 0.4Af, h/2 and h/A and refer to the nonrotating mode 
A of Table 1 when evolved in a fixed spacetime (left panel) and in a dynamical one (right 
panel). Indicated with a dot-dashed line is the slope for a second-order convergence rate. This 
helps to see that both formulations yield an almost second-order convergent measure of the 
eigenfrequencies of the F-mode, with the new formulation having a truncation error which 
is always smaller than the one coming from the standard formulation. Given the importance 
of an accurate measurement of the eigenfrequencies to study the mode properties of compact 
stars, we believe that Figs. 7 and 8 provide an additional evidence of the advantages of the 
new formulation. 

Finally, we note that a behaviour similar to the one shown in Fig. 7- 8 has been found also for 
rotating stars although in this case the comparison is possible only for evolutions within the 
Cowling approximation since we lack a precise perturbative estimate of the eigenfrequency 
for model B of Table 1 for a dynamical spacetime. 



4.4. Cylindrical Shock Reflection 

One of the most important properties of HRSC schemes is their capability of handling the 
formation of discontinuities, such as shocks, which are often present and play an important 
role in many astrophysical scenarios. Tests involving shocks formation are usually quite 
demanding and codes that are not flux-conservative can also show numerical instabilities 
or difficulties in converging to the exact solution of the problem. Since both the new and 
the standard formulation solve the relativistic hydrodynamics equations as written in a flux- 
conservative form, they are both expected to be able to correctly resolve the formation of 
shocks, although each with its own truncation error. In the following test we consider one 
of such discontinuous flows and show that the new formulation provides a higher accuracy 
with respect to the standard one, stressing once again the importance of the definition of the 
conserved variables. 

More specifically, we consider a one-dimensional test, first proposed by [51], describing the 
reflection of a shock wave in cylindrical coordinates. The initial data consist of a pressureless 
gas with uniform density po = 1-0, radial velocity = 0.999898, corresponding to an initial 
Lorentz factor Wq = 70.0 and an internal energy which is taken to be small and proportional 
to the initial Lorentz factor, i.e., e = 10~^(Wo). During the evolution an ideal-fluid EOS (30) 
is used with a fixed adiabatic index F = 4/3. The symmetry condition at x = produces 
a compression and generates an outgoing shock in the radial direction. The analytic solution 
for the values of pressure, density, gas and shock velocities are given in [51]. From them one 
can determine the position of the shock front at any time t 

^ {T-l)W,H\ 

Wa + l ^ ' 

This can then be used to compare the accuracy of the two formulations. 
In the left panel of Fig. 9 we show the value of the radial component of the velocity as 
a function of x at a time t = 0.002262 ms and for a resolution of Ii/Mq = 6.25 x 10~^. 
The solid line represent the analytic solution, the short-dashed line the numerical solution 
computed with the new formulation and the long-dashed line the one obtained with the 
standard formulation. As it is evident from the inset, the position of the shock is very well 
captured by both formulations, but the new one is closer to the exact one at this time. 
To compare with the exact prediction given by expression (33), we compute the numerical 
position of the shock as the middle of the region where the value of the velocity moves 
from the pre-shock value w+ to the post-shock one v~ (in practice, we fit a straight line 
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Figure 9. Left panel: Comparison of the velocity profiles for the two formulations in the 
solution of the axisymmetric shock-tube test with a resolution ofh = 6.25 X 10^^ A/q The 
solid line shows the exact position after a time f = 0.002262 ms, while the short-dashed 
and the long-dashed lines represent the solutions with the new and the standard formulations, 
respectively. Right panel: Comparison of the error in the determination of the position of 
the shock in the two formulations. Note the first-order convergence rate as expected for 
discontinuous flows. 



between the last point of the constant post-shock state and the first point of the constant 
pre-shock state and evaluate the position at which this function has value (w+ + v^)/2.). 
The right panel of Fig. 9, shows the relative error 1 — (xg)num/(2;s)anai in the position of 
the shock at time i = 0.002622 ms and for five different resolutions: h = 0.01 Mq, h/8, 
h/AO, h/80 and /i/160. Indicated with a dashed line is the error computed when using the 
standard formulation, while indicated with a solid line is the error coming from the new 
formulation. Note that both formulations show a first-order convergence, as expected for 
HRSC schemes in the presence of a discontinuous flow, but, as for the other tests, also in this 
case the new formulation has a smaller truncation error. A similar behaviour is shown also by 
other quantities in this test but these are not reported here. 

It is useful to note that the difference between the two formulations in this test is smaller than 
in the previous ones, being of a factor of a few only and not of orders of magnitude. We 
believe this is due in great part to the fact that, in contrast with what happens for stars, the 
solution in the most troublesome part of the numerical domain {i.e. near a: ~ 0, z ~ 0) is not 
characterized by particularly large values of the fields or of the fluid variables. In support of 
this conjecture is the evidence that at earlier times, when the shock is closer to the axis, both 
the absolute errors and the difference between the two formulations are larger 

5. Conclusion 

A number of astrophysical scenarios can be very conveniently studied numerically by 
assuming they possess and preserve a rotation symmetry around an axis. Such an assumption 
reduces the number of spatial dimensions to be considered and thus the computational costs. 
This, in turn, allows for a more sophisticated treatment of the physical and astrophysical 
processes taking place and, as a result, for more realistic simulations. 
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We have presented a new numerical code developed to solve in Cartesian coordinates the 
full set of general relativistic hydrodynamics equations in a dynamical spacetime and in 
axisymmetry. More specifically, the new code solves the Einstein equations by using the 
"cartoon" method, while HRSC schemes are used to solve the hydrodynamic equations 
written in a conservative form. An important feature of the code is the use of a novel 
formulation of the equations of relativistic hydrodynamics in cylindrical coordinates. More 
specifically, by exploiting a suitable definition of the conserved variables, we removed 
from the source of the flux-conservative equations those terms that presented coordinate 
singularities at the axis and that are usually retained in the standard formulation of the 
equations. Despite their simplicity, the changes made to the standard formulation can produce 
significant improvements on the overall accuracy of the simulations with a truncation error 
which is often several orders of magnitude smaller 

In order to assess the validity of the new formulation and compare its accuracy with that of 
the formulation which is most commonly used in Cartesian coordinates, we have performed 
several tests involving the evolution of oscillating spherical and rotating stars, as well as 
shock-tube tests. In all cases considered we have shown that the codes implementing the 
two formulations yield the expected convergence rate but also that the new formulation is 
always more accurate, often considerably more accurate, than the standard one. 
In view of its simplicity, the new formulation of the equations can be implemented 
straightforwardly in codes written using the standard formulation and we recommend its use 
for all situations in which an axisymmetric problem needs to be investigated in full general 
relativity and in Cartesian coordinates. 
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6. Appendix A 

In what follows we recall the essential features of the "cartoon " method for the solution of the 
field equations in Cartesian coordinates. Consider therefore the computational domain to have 
extents < x,z < dmax and — Aj/ < y < Ay, where dmax refers to the location of the outer 
boundary. Reflection symmetry with respect to the z = 0-plane can additionally be assumed. 
The Einstein equations are then solved only on the y ~ 0-plane with the derivatives in the y- 
direction being computed with second-order centred stencils using the points at — Ay, , Ay. 
Taking into account axisymmetry, the rotation in the (.t, y) plane is defined as 



where R{(f>)^^ = R{—4>) and the rotation angle is defined as,(j) — tan^^ (±Aj// yjx^ + (Ay)^). 
As commented in the main text, the values of all the quantities on the ±Ay planes are 
computed via rotations of the corresponding values on the y = 0-planes. More specifically. 




- sin(0) 
cos(0) 
1 




(34) 
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the components of an arbitrary vector field on the ±Ay planes are computed via a (f)- 
rotation as 

= tJo) cos{4>) - ) s\n{4>) , (35) 
Ty = Ti") sin(0) + T^") cos(0) , (36) 
T, = ri") , (37) 

where Tj'^-* denote the corresponding components at {y/ + (A?/)^, 0, z), which are 
computed via a Lagrangian interpolation from the neighboring points on the z-axis. Similarly, 
the components of an arbitrary tensor field tensor will be computed as 

T,, = Tf} cos^ (0) - Ti°) sin(2^) + Tg^ sin^ (0) , (38) 

T.y = sin(20) - cos(2(/.) + ^T^^ sin(20) , (39) 

T,, = T^) siii^ (0) - Ti;;) sin(2(/.) + T^^j) cos^ (</.) , (40) 

T,, Ti^) cos(0) ~ 7^(0) sin(0) , (41) 

Ty, = Ti^) sin(0) + Ty^o) cos((/)) , (42) 

T,, = ri°) . (43) 
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